
#include "glodef.h"
#include "errlog.h"
#include "fe.h"
#include "fe_matr.h"

#include <stdlib.h>
#include <stdio.h>

#define POISSON_ID           0x3F1EC68B

typedef struct
{
    unsigned id;
    unsigned eltype;
    int null_dof;
    grid3d_t *grid;
    fe_matr_p poisson;
} poisson_t;

static poisson_t *get_pdata(void *data)
{
    poisson_t *p=data;
    if((data==NULL)||(p->id!=POISSON_ID))
        error_msg("poisson.c: invalid poisson_t structure passed.\n");
    return p;
}




static int null_dof= -1;

static void eo_poisson(const grid3d_t *grid, unsigned eltype, const int *dofs, unsigned **eos)
{
    static unsigned _eo[MAX_LDOFS];
    int nn, i;
    if(eos)
        *eos=_eo;
    else
        return;
    if(null_dof== -1)
        error_msg("technical error 11 (0x0B).\n");
    nn=fe_dofs(grid, 0, eltype, NULL);
    for(i=0; i<nn; i++)
        _eo[i]=dofs[i]==null_dof?eo_bc:eo_add;
}

int SolvePoisson(void *pdata, scalar U, scalar F, scalar T,
                  unsigned F_tested, unsigned reassemble)
{
    poisson_t *pd=get_pdata(pdata);
    if(!F_tested)
    {
        test_scalar(pd->grid, pd->eltype, F, T);
        copy_scal(T, F);
    }
    set_value_scal(F, pd->null_dof, 0.0);
    set_value_scal(U, pd->null_dof, 0.0);

    assemble_scal(F);
    assemble_scal(U);

    if(reassemble)
    {
        null_dof=pd->null_dof;
        fe_matr_assemble(pd->grid, 1, &pd->poisson);
        null_dof= -1;
    }
    PrintLog("Solving Poisson problem.\n");
    int ret = fe_matr_solve_scalar(pd->poisson, U, F, T, 0.0);
    if(!ret)
        warning_msg("Poisson solver did not converge.\n");
    return ret;
}

void ReassemblePoisson(void *pdata)
{
    poisson_t *pd=get_pdata(pdata);
    fe_matr_assemble(pd->grid, 1, &pd->poisson);
}

void *InitPoisson(grid3d_t *grid, unsigned eltype)
{
    LOG_MPI_WTIME_START;

    poisson_t *pd;
    InitFE(1, &eltype, grid);

    PrintLog("Initialize Poisson solver (%s) begin.\n", __func__);

    pd=malloc(sizeof(poisson_t));
    if(pd==NULL)
        error_msg("poisson.c: memory allocation error.\n");
    pd->id=POISSON_ID;
    pd->eltype=eltype;
    {
        int ndof;
        ndof=fe_ndof(grid, eltype);
        for(pd->null_dof=0; pd->null_dof<ndof; pd->null_dof++)
            if(fe_dof_bcond(grid, eltype, pd->null_dof)>BC_NON)
                break;
    }
    null_dof=pd->null_dof;

    pd->grid=grid;
    pd->poisson=fe_matr_alloc(eltype, local_stif, eo_poisson, grid, TRUE,
#if defined(USE_CSC)  
            "ss"
#else
            NULL
#endif
            );
    fe_matr_assemble(pd->grid, 1, &pd->poisson);
    null_dof= -1;

    /*
    fe_matr_t *m=pd->poisson;
    static char fname[100];
    extern int curr_ti;
    sprintf(fname, "matr/poisson_%03d.txt", curr_ti);
    PrintSMatrix(m->matr, fname);
     */

    LOG_MPI_WTIME_STOP;

    return pd;
}

void DonePoisson(void *pdata)
{
    poisson_t *pd=get_pdata(pdata);
    fe_matr_free(pd->poisson);
    free(pd);
    PrintLog("Poisson solver done.\n");
}


